%%LOOPSTACK%%
%Iterates through variables for acoustics=>articulation reconstruction.
%These variables include:
%Lag: the difference in real time between ema and acoustic samples
%considered to be at the same 'time'
%Stagger: The distance in samples between central window and two side
%windows of acoustic features concatenated into one feature vector
%Central window size
%Side window size
%number of principle components

%Requires an m x n array of acoustic features called ACOUSTICS where each 
%row represents a file and each column a window size. Each cell contains a
% t x u matrix of u acoustic features at t times. (see
% feature_extraction.m)

%Requires an m x n array of EMA features called EMA where each row
%represents a file and each column an EMA feature space. Each cell contains
%a t x u matrix of u ema features at t times

%Requires an m x 1 array of annotations (see Praat Script texttablewriter)
%for m files. Each cell contains a t x 2 matrix of t utterances' time of
%onset (col 1) and offset (col 2)


%declaring range of lag and stagger variables in samples (200 Hz)

LAG = linspace(-20, 20, 41);
STAGGER = linspace(1, 40, 40);
rotation = 5;
ema_features_raw = 15;
ema_features_palate = 9;
ema_features_td = [];

ema_features = [ema_features_raw, ema_features_palate, ema_features_td];



METHODS = {'t = multlr(ema_training, acoustics_training(:, 1:n)); ema_result = transp(transp(t)*transp(acoustics_test(:, 1:n)));',...
           [ ...
           'net = mlp(1, 3, 1, ''linear'');'...
           'options = zeros(1,18); options(1) = 0;options(14) = 100;'...
           '[net, options] = netopt(net, options, acoustics_training(:, 1:n), ema_training(:, 1:n), ''scg'');   ']};
       
EXP_PARAMS = {{''},...
              {'nin = 60; nhidden = 100; nout = 9; alpha = 0.01;',...
              'nin=60; nhidden=150; nout=9; alpha = 0.01'}}; 

%%Add in other methods (here is MLR) as commands as strings. 
%%They are evaluated in a loop at the bottom
%At minimum, method must output a m x n matrix of n reconstructed
%ema features at m time samples called 'ema_result';


%the results need to be in separate matrices since different ema
%feature spaces have a different number of features

ema_raw = struct('COR', [], 'ERROR', []);
ema_palate = struct('COR', [], 'ERROR', []);
ema_tada = struct('COR', [], 'ERROR', []);

acoustic_windows = length(ACOUSTICS);
[n_dataset_samples, acoustic_features] = size(ACOUSTICS{1});
acoustic_features = acoustic_features*3;  %% multiplied by 3 for the two staggered samples
randomizing_vector = rand([n_dataset_samples, 1]);
    
SEGS = unique(labels);



%%Starting down the loop

for g = 1:3
    
    COR = cell(numel(LAG), numel(STAGGER), acoustic_windows, acoustic_windows, rotation, acoustic_features, length(METHODS), ema_features(g));  
    ERR = cell(numel(LAG), numel(STAGGER), acoustic_windows, acoustic_windows, rotation, acoustic_features, length(METHODS), ema_features(g));
    
    for h = 1:length(SEGS)
        
        %gets the indexes for the hth segment
        seg_i = labels_i((find((strcmp(labels, SEGS{h}))==1)));
        n_samples = length(seg_i);
        n_samples = (floor(n_samples/rotation))*rotation; %this just ensures that
        seg_i = seg_i(1:n_samples);
        n_test_samples = n_samples/rotation;
        n_training_samples = n_samples-n_test_samples;
        

        %assmbling matrix of ema features
        ema = EMA{g}(seg_i, :);
    
        for i = 1:numel(LAG)
            lag = LAG(i);
            for j = 1:numel(STAGGER)
                stagger = STAGGER(j);
                for k = 1:acoustic_windows
                    for l = 1:acoustic_windows
                        acoustics = horzcat(ACOUSTICS{k}(seg_i+lag, :), ACOUSTICS{l}(seg_i+lag+stagger, :), ACOUSTICS{l}(seg_i+lag-stagger, :));
                        acoustics = acoustics(1:n_samples, :);
                        
                        %principle components analysis
                        coefs = princomp(acoustics);
                        scores = (coefs*acoustics')';
                    
                        % randomizing the acoustics & ema
                       
                        acoustics = horzcat(scores, randomizing_vector(seg_i));
                        acoustics = sortrows(acoustics, acoustic_features+1);
                        acoustics = acoustics(:, 1:acoustic_features);
                    
                        ema = horzcat(ema, randomizing_vector(seg_i));
                        ema = sortrows(ema, ema_features(g)+1);
                        ema = ema(:, 1:ema_features(g));
               
                        for m = 1:rotation
                            
                            %assembling test/training matrices
                            
                            test_i = ((m-1)*n_test_samples)+1:m*n_test_samples;
                            training_i = setdiff(1:n_samples, test_i);
                            acoustics_training = acoustics(training_i, :);
                            acoustics_test = acoustics(test_i, :);
                            ema_training = ema(training_i, :);
                            ema_test = ema(test_i, :);
                        
                            for n = 1:acoustic_features
                                
                                for o = 1:length(EXP_PARAMS)  % METHOD index is implicit
                                    for e = 1:length(EXP_PARAMS{o})

                                        eval(EXP_PARAMS{o}{e});   % TODO
                                        eval(METHODS{o});


                                        for p = 1:ema_features(h)
                                           COR{h, i, j, k, l, m, n, o, p} = corr(ema_result(:, p), ema_test(:, p)); 
                                           ERR{h, i, j, k, l, m, n, o, p} = mean((ema_result(:, p) - ema_test(:, p)).^2);
                                            
                                            %h = segment
                                            %i = lag
                                            %j = stagger
                                            %k = central window size
                                            %l = side window size
                                            %m = rotation
                                            %n = number of principle components
                                            %o = method
                                            %p = ema features
                                            
                                        end
                                    end
                                end
                            end 
                        end
                    end
                end
            end
        end
    end
    
    if g == 1
        ema_raw.COR = COR;
        ema_raw.ERR = ERR;
    elseif g == 2
        ema_palate.COR = COR;
        ema_palate.ERR = ERR;
    else 
        ema_tada.COR = COR;
        ema_tada.ERR = ERR;
    end
    
end